###iso
setwd("D:/入侵/数据/分析2021/6 iso")
library(SIBER)
library(dplyr)
library("reshape2")
iso_dat<-read.csv("iso_CN.csv",header=TRUE)
iso_C<- iso_dat %>% 
  filter(iso == "C") %>% 
  melt(.,  measure.vars = c("A", "B", "C","D","E"), 
       variable.name='Site', value.name='C_val', na.rm = TRUE)
iso_N<- iso_dat%>% 
  filter(iso == "N") %>% 
  melt(., measure.vars = c("A", "B", "C","D","E"),
       variable.name='Site', value.name='N_val', na.rm = TRUE)
mydata=data.frame(iso1=iso_C$C_val,iso2=iso_N$N_val,group=iso_C$func,community=iso_C$Community)
mydata$community = factor(mydata$community, levels=c('O','R','P','S','PS'))

mydata=na.omit(mydata)  #Delete missing values
#mydata=read.csv("consumers2.csv",header=TRUE)
siber.example <- createSiberObject(mydata)  

#Calculate sumamry statistics for each group: TA, SEA and SEAc
group.ML <- groupMetricsML(siber.example) 
write.csv(group.ML,file = "group.ML.csv",row.names = T)

community.ML <- communityMetricsML(siber.example)
write.csv(community.ML,file = "community.ML.csv",row.names = T)

library(rjags)
parms <- list()
parms$n.iter <- 2 * 10^4   # number of iterations to run the model for
parms$n.burnin <- 1 * 10^3 # discard the first set of values
parms$n.thin <- 10     # thin the posterior by this many
parms$n.chains <- 2        # run this many chains

priors <- list()
priors$R <- 1 * diag(2)
priors$k <- 2
priors$tau.mu <- 1.0E-3

ellipses.posterior <- siberMVN(siber.example, parms, priors)

SEA.B <- siberEllipses(ellipses.posterior)
siberDensityPlot(SEA.B, xticklabels = colnames(group.ML),
                 xlab = c("Community | Group"),
                 ylab = expression("Standard Ellipse Area " ('\u2030' ^2) ),
                 ylims = c(0,100),
                 bty = "L",
                 las = 1,
                 main = "SIBER ellipses on each group"
) #28*5

cr.p <- c(0.95, 0.99) 

SEA.B.credibles <- lapply(
  as.data.frame(SEA.B),
  function(x,...){tmp<-hdrcde::hdr(x)$hdr},
  prob = cr.p)
write.csv(SEA.B.credibles,file = "SEA.B.credibles.csv",row.names = T)

SEA.B.modes <- lapply(
  as.data.frame(SEA.B),
  function(x,...){tmp<-hdrcde::hdr(x)$mode},
  prob = cr.p, all.modes=T)
points(1:ncol(SEA.B), SEA.B.modes, col="red", pch = "x", lwd = 2)
write.csv(SEA.B.modes,file = "SEA.B.modes.csv",row.names = T)

mu.post <- extractPosteriorMeans(siber.example, ellipses.posterior)

layman.B <- bayesianLayman(mu.post)

layman_list <-  list(1,2,3,4,5)
for (i in layman_list) {
  pdf(paste('layman_', i, '.pdf', sep = ''), width = 5, height = 5)
  siberDensityPlot(layman.B[[i]], xticklabels = colnames(layman.B[[i]]),
                   bty="L", ylim = c(0,13))
  dev.off()
  credibles <- lapply(
    as.data.frame(layman.B[[i]]),
    function(x,...){tmp<-hdrcde::hdr(x)$hdr},
    prob = cr.p)
  write.csv(credibles,file = paste("credibles_", i, '.csv',sep = ""),row.names = T)
  
  modes <- lapply(
    as.data.frame(layman.B[[i]]),
    function(x,...){tmp<-hdrcde::hdr(x)$mode},
    prob = cr.p, all.modes=T)
  write.csv(modes,file = paste("modes_", i, '.csv',sep = ""),row.names = T)
}


##Figure drawing
library(MixSIAR)
library(dplyr)
library("reshape2")
iso_dat<-read.csv("iso_CN.csv",header=TRUE)
com = list("O","R","P","S","PS")
for (i in com) {
  iso_i=iso_dat %>% filter(Community == i)
  iso_C<- iso_i %>% 
    filter(iso == "C") %>% 
    melt(.,  measure.vars = c("A", "B", "C","D","E"), 
         variable.name='Site', value.name='C_val', na.rm = TRUE)
  iso_N<- iso_i%>% 
    filter(iso == "N") %>% 
    melt(., measure.vars = c("A", "B", "C","D","E"),
         variable.name='Site', value.name='N_val', na.rm = TRUE)
  mydata=data.frame(Group=iso_C$func,d13C=iso_C$C_val,d15N=iso_N$N_val,community=iso_C$Community)
  write.csv(mydata,file = paste("com_", i, '.csv',sep = ""),row.names = F)
  
  mydata=read.csv(paste("com_", i, '.csv',sep = ""),header=T)
  mydata=data.frame(iso1=mydata$d13C,iso2=mydata$d15N,group=mydata$Group,community=mydata$community)
  siber.example <- createSiberObject(mydata)
  
  dat_5=mydata%>% filter(group=="5")
  dat_6=mydata%>% filter(group=="6")
  dat_7=mydata%>% filter(group=="7")
  dat_8=mydata%>% filter(group=="8")
  dat_P=mydata%>% filter(group=="P")
  dat_S=mydata%>% filter(group=="S")
  #tmp<-iso_data2 %>% filter(Community == i,func=="P")
  
  plot(x=dat_5$iso1, y=dat_5$iso2, pch = 20,col="#F06C79",
       xlab = expression({delta}^13*C~'\u2030'), ylab = expression({delta}^15*N~'\u2030'),
       xlim = c(-35,-10),ylim = c(0,25))
  points(x=dat_6$iso1, y=dat_6$iso2, pch = 20,col="#4DB34D")
  points(x=dat_7$iso1, y=dat_7$iso2, pch = 20,col="#FDCF45")
  points(x=dat_8$iso1, y=dat_8$iso2, pch = 20,col="#336FCC")
  points(x=dat_P$iso1, y=dat_P$iso2, pch = 20,col="#000000")
  points(x=dat_S$iso1, y=dat_S$iso2, pch = 20,col="#000000")
  
  plotGroupEllipses(siber.example, n = 100, p.interval = 0.95, ci.mean = T, lty = 1, lwd = 2)
}
#F06C79   5 red--black
#4DB34D   6 green--red
#FDCF45   7 yellow--green
#336FCC   8 blue--blue
##          black--brilliant blue


for (i in com) {
  iso_i=iso_dat %>% filter(Community == i)
  iso_C<- iso_i %>% 
    filter(iso == "C") %>% 
    melt(.,  measure.vars = c("A", "B", "C","D","E"), 
         variable.name='Site', value.name='C_val', na.rm = TRUE)
  iso_N<- iso_i%>% 
    filter(iso == "N") %>% 
    melt(., measure.vars = c("A", "B", "C","D","E"),
         variable.name='Site', value.name='N_val', na.rm = TRUE)
  mydata=data.frame(Group=iso_C$func,d13C=iso_C$C_val,d15N=iso_N$N_val,community=iso_C$Community)
  write.csv(mydata,file = paste("com_", i, '.csv',sep = ""),row.names = F)
  
  mix.filename <- paste("com_", i, '.csv',sep = "")
  mix <- load_mix_data(filename=mix.filename,
                       iso_names=c("d13C","d15N"),
                       factors="Group",
                       fac_random=FALSE,
                       fac_nested=FALSE,
                       cont_effects=NULL)

  source <- load_source_data(filename="sources.csv",
                             source_factors=NULL,
                             conc_dep=FALSE,
                             data_type="raw",
                             mix)
  
  discr <- load_discr_data(filename="discr.csv", mix)
  
  plot_data(filename=paste("isospace_plot_", i, sep = ""),
            plot_save_pdf=TRUE,
            plot_save_png=FALSE,
            mix,source,discr)
  
  model_filename <- paste("MixSIAR_model_", i, '.txt',sep = "")
  resid_err <- TRUE
  process_err <- FALSE
  write_JAGS_model(model_filename, resid_err, process_err, mix, source)
  
  #jags.1 <- run_model(run="test", mix, source, discr, model_filename, alpha.prior=1, resid_err, process_err)
  jags.1 <- run_model(run="normal", mix, source, discr, model_filename, alpha.prior=1, resid_err, process_err)
  
  output_options <- list(summary_save = TRUE,                 
                         summary_name = paste("summary_statistics_", i, sep = ""), 
                         sup_post = FALSE,                    
                         plot_post_save_pdf = FALSE,           
                         plot_post_name = "posterior_density",
                         sup_pairs = FALSE,             
                         plot_pairs_save_pdf = FALSE,    
                         plot_pairs_name = "pairs_plot",
                         sup_xy = FALSE,           
                         plot_xy_save_pdf = FALSE,
                         plot_xy_name = "xy_plot",
                         gelman = FALSE,
                         heidel = FALSE,  
                         geweke = FALSE,   
                         diag_save = FALSE,
                         diag_name = "diagnostics",
                         indiv_effect = FALSE,       
                         plot_post_save_png = FALSE, 
                         plot_pairs_save_png = FALSE,
                         plot_xy_save_png = FALSE)
  output_JAGS(jags.1, mix, source, output_options)
  graphics.off()
}